Uncovering the phenotypic consequences of multi-locus imprinting disturbances using genome-wide methylation analysis in genomic imprinting disorders

Imprinted genes are regulated by DNA methylation of imprinted differentially methylated regions (iDMRs). An increasing number of patients with congenital imprinting disorders (IDs) exhibit aberrant methylation at multiple imprinted loci, multi-locus imprinting disturbance (MLID). We examined MLID and its possible impact on clinical features in patients with IDs. Genome-wide DNA methylation analysis (GWMA) using blood leukocyte DNA was performed on 13 patients with Beckwith–Wiedemann syndrome (BWS), two patients with Silver–Russell syndrome (SRS), and four controls. HumanMethylation850 BeadChip analysis for 77 iDMRs (809 CpG sites) identified three patients with BWS and one patient with SRS showing additional hypomethylation, other than the disease-related iDMRs, suggestive of MLID. Two regions were aberrantly methylated in at least two patients with BWS showing MLID: PPIEL locus (chromosome 1: 39559298 to 39559744), and FAM50B locus (chromosome 6: 3849096 to 3849469). All patients with BWS- and SRS-MLID did not show any other clinical characteristics associated with additional involved iDMRs. Exome analysis in three patients with BWS who exhibited multiple hypomethylation did not identify any causative variant related to MLID. This study indicates that a genome-wide approach can unravel MLID in patients with an apparently isolated ID. Patients with MLID showed only clinical features related to the original IDs. Long-term follow-up studies in larger cohorts are warranted to evaluate any possible phenotypic consequences of other disturbed imprinted loci.


Introduction
Genomic imprinting is an epigenetic regulatory mechanism, leading to a parent-of-origin-specific expression of a small subset of genes [1]. The life cycle of the genomic imprints in mammals consists of three stages: establishment of parental imprinting marks in the germline during gametogenesis, imprint maintenance through fertilization and early embryonic Amsterdam, Netherlands) following the manufacturer's instructions. Targeted bisulfite pyrosequencing assays covering four and seven consecutive CpG sites for ICR1 and ICR2, respectively, were performed using a PyroMark Q24 pyrosequencer (Qiagen, Hilden, Germany) (S1 Table). Methylation levels were calculated as a percentage of methylated cytosine [% mC = mC/(mC+C)] for each CpG site using PyroMark Q24 Software (v.1.0.10; Qiagen). The standard deviation scores (SDS) of altered DNA methylation level (%mC) at ICR1 or ICR2 in patients with BWS and SRS was calculated based on the average DNA methylation level detected in 20 samples from age-and sex-matched controls.
For patients with BWS, the clinical score was calculated using the BWS consensus scoring system: 2 points per each cardinal features (macroglossia, exomphalos, lateralized overgrowth, multifocal Wilms tumor, prolonged hyperinsulinism, distinct pathologic findings), and 1 point per each suggestive features (large for gestational age, facial naevus flammeus, polyhydramnios or placentomegaly, ear creases or pits, transient hypoglycemia, embryonal tumors, nephromegaly or hepatomegaly and umbilical hernias or diastasis recti) [6]. Patients with score of �4 were regarded as classical BWS. For patients with SRS, the clinical score was calculated using the Netchine-Harbison clinical scoring system (NH-CSS), including small for gestational age, postnatal growth failure, relative macrocephaly at birth, protruding forehead, body asymmetry, feeding difficulties and/or low body mass index [7]. Age-and sex-specific SDS for serial growth parameters including height, weight, head circumference, and/or body mass index were assigned based on Fenton growth references at birth [22] and 2017 Korean National Growth Charts at postnatally [23]. Overgrowth or undergrowth were defined as height or weight greater than or less than two standard deviations compared to age-and sexmatched controls [6,7]. Developmental and pubertal status was evaluated by a clinical geneticist and an endocrinologist, respectively. Skeletal abnormalities, including leg length discrepancy, scoliosis, and thoracic deformity, were assessed by an orthopedic surgeon using x-ray images. Clinical data for other phenotypes, including cardiac, gastrointestinal, and central nervous system anomalies, were retrospectively collected from medical records.

GWMA and data analysis
A peripheral blood sample (5 mL) was collected in each participant at initial visit, and genomic DNA was extracted from all types of peripheral blood leukocytes including granular and agranular leukocytes using the DNA isolation kit (QIAGEN, Hilden, Germany). GWMA was performed in 13 patients with BWS, two patients with SRS, and four healthy children controls. Sodium bisulfite conversion of genomic DNA was conducted using an EZ DNA Methylation Kit (Zymo Research). After bisulfite conversion, each sample was whole-genome amplified, enzymatically fragmented, and applied to the Illumina HumanMethylation850 BeadChip (Illumina, Inc., San Diego, CA, USA), which contains probes at more than 850,000 CpG loci. The amplified DNA was annealed to allele-specific primers linked to two individual bead types, which correspond to each CpG locus-one methylated (C) and one unmethylated (T). After single-base extension using DNP-and biotin-labeled ddNTPs, the array was fluorescently labeled, followed by washing and coating. The Illumina BeadArray Reader, and Illumina Gen-omeStudio v. 2011.1 were used for image reading and extracting image intensities. The location of the CpG relative to UCSC CGI was assigned to each probe; shores (sequences of 0-2 kb from the CGI), shelves (sequences 2-4 kb from the CGI), and open sea (sequences located outside these regions) [24]. The CpGs were functionally classified as proximal promoters (CpGs located within 200 or 1500 bp upstream of transcription start sites, exon 1, and 5 0 UTRs), 3 0 UTRs, gene bodies, and intergenic regions [24]. Arrays were processed at Macrogen (Seoul, Korea) by following the manufacturer's instructions.
The β values were extracted for each CpG by subtracting background intensity computed from negative controls and taking the ratio of the methylated signal intensity to the combined intensity of methylated and unmethylated signals ranging from 0 (no methylation) to 1 (100% methylation). From 865,918 CpGs in total, 865,544 CpGs were detected by detection P-value of <0.05. As a quality control step, we eliminated probes associated with single-nucleotide polymorphism, cross-hybridization probes, and probes that corresponded to the X or Y chromosome from the dataset, leaving 745,915 CpGs to be analyzed. The Beta Mixture Quantile (BMIQ) normalization of filtered data was conducted to reduce assay bias [25]. To confirm the abnormal methylation levels of disease-specific iDMRs and detect MLID in each patient, the Crawford-Howell t-test was performed using preprocessed data [26]. For 809 CpGs on 77 iDMRs defined by Monk [27] and Joshi [28], the DNA methylation level difference (delta beta, Δβ) was calculated by subtracting mean β value of the control group from β value of each patient. Differentially methylated position (DMP) was defined as the absolute value of Δβ (| Δβ|) of >0.2 (>20% difference in DNA methylation) and P-value of <0.05. Aberrantly methylated iDMRs were defined as iDMRs showing differentially methylated levels in at least two consecutive probes within iDMRs which included at least four significant probes [19]. The VTRNA2-1:DMR was considered as insignificant because it has been shown to be polymorphically imprinted in general population, showing frequent hypomethylation with various degrees [29]. We considered each patient as having MLID if showing aberrant methylation in one or more additional iDMRs other than the disease-related iDMRs.

Exome sequencing
Three samples from patients with BWS exhibiting MLID were subjected to whole-exome sequencing (Fig 1). The exome sequencing data of the patient were analyzed using SureSe-lectXT Human All Exon V5 (Agilent Technologies, CA, USA), and the library was prepared according to the manufacturer's protocol. An Illumina HiSeq2500 platform was used for library sequencing with 2 × 151 paired-end reads. Reads were aligned to the human genome build 37 (hg19) using Burrows-Wheeler Aligner (v. 0.7.17), and PCR duplicates were filtered out using Picard software (v. 2.9.0). The binary alignment map (BAM) file was realigned and recalibrated using SAMtools (v. 1.9), and the Genome Analysis Toolkit (v. 4.1.2). Variants were called using HaplotypeCaller and annotated using SnpEff, ANNOVAR, and InterVar. Sequence variants were compared with databases such as Genome Aggregation Database (gno-mAD), Human Gene Mutation Database, and ClinVar. For further analysis, variants with zero allele frequency in the gnomAD were selected for autosomal dominant genes, and variants with allele frequency <0.01% were selected for autosomal recessive genes. The low-frequency variants (allele frequency 0.05-0.25) were detected using MuTect2. Copy number variation and low-frequency variants were analyzed as described previously [30]. After variant classification using guidelines of the American College of Medical Genetics and Genomics guidelines, we considered pathogenic or likely pathogenic variants as causative variants [31]. Variant screening for candidate genes was performed and genes involved in imprinting establishment and maintenance during embryonic development were analyzed (S2 Table) [32,33].

Statistical analysis
Descriptive data are presented as mean ± standard deviation, and categorical variables as counts and proportions. The group differences according to the MLID (mono-locus vs. multilocus) were investigated using the Mann-Whitney U test and Fisher's exact test. Statistical analysis was performed using SPSS v. 25.0 (SPSS Inc. Chicago, IL, USA) and R software v. 3.6.0 (The Comprehensive R Archive Network, Vienna, Austria; https://cran.r-project.org). A P-value of <0.05 was considered as statistically significant. Table 1 shows the clinical and molecular characteristics of 13 patients with BWS (six males, seven females), and two patients with SRS (one male, one female). The mean follow-up duration was 4.6 ± 3.2 years in BWS, and 9.0 ± 2.3 years in SRS. Patients with BWS comprised two different epigenotypes (12 IC2-LoM, one IC1-GoM) and all met the criteria for classical BWS (a score of �4). The mean SDS for height and weight at initial visit (at a mean age of 1.8 ± 2.4 years) in patients with BWS were +2.2 ± 1.6, and +1.9 ± 1.1, respectively. Two patients with SRS had the IC1-LoM epigenotype and both met more than four out of six NH-CSS features. The SDS for height and weight at initial visit were -4.2, and -3.4 in SRS1 (at 3.6 months of age), and -4.1, and -7.2 in SRS2 (at 4.1 years of age), respectively. Both patients with SRS were treated with growth hormone, starting at a mean age of 4.3 years at a dosage for small for gestational age (0.3 mg/kg/week). SRS2 received combination therapy with a gonadotropin-releasing hormone agonist due to advanced puberty.

Genome-wide alteration of DNA methylation
Hierarchical cluster analysis of the methylated region using GWMA identified 74 significant DMPs in 15 patients with ID compared with four healthy children controls: 29 (39.2%) hyper DMPs and 45 (60.8%) hypo DMPs. Patients with ID were generally classified into two categories based on the genome-wide methylation pattern: (1) category 1 of 12 patients with BWS (all with IC2-LoM); (2) category 2 of two patients with SRS (both with IC1-LoM) and one patient with BWS (with IC1-GoM) (Fig 2). In a single-sample analysis for each patient, the number of significantly hypermethylated DMPs per patient ranged from 175 (BWS8) to 17,035 (BWS7), while the number of significantly hypomethylated DMPs ranged from 375 (SRS2) to 3,760 (BWS7) ( Table 2).
When focusing on 809 CpGs in 77 iDMRs defined by Monk et al. [27] and Joshi et al. [28], a total of 232 CpGs in 17 iDMRs (two paternally and 15 maternally methylated) were identified as differentially methylated in patients with ID, compared with the control group (mean methylation level of 0.4-0.6). All patients with BWS and SRS showed aberrant methylation in disease-related iDMRs, which were consistent with those identified in targeting analysis using bisulfite pyrosequencing (Table 1). In addition, the GWMA data revealed 14 additional hypomethylation other than the disease-related iDMRs, suggesting MLID in four patients: three patients with BWS (BWS8, BWS9, and BWS12), and one patient with SRS (SRS1) (S3 Table).

Clinical features of patients with and without MLID
All patients with MLID showed the classical phenotype of BWS or SRS without any clinically distinguishable features ( technique), parental age at conception, prenatal and postnatal growth parameters, BWSrelated phenotypes, and clinical score (Table 3).

Exome analysis of putative MLID-causative variants
In order to investigate the putative MLID-causative variants, whole-exome sequencing using blood leukocyte DNA was performed in three patients with BWS who exhibited multiple  hypomethylation. No copy number variant was found in the exome data of three patients. In addition, no causative variants were identified among candidate genes, including KHDC3L, NLRP2, NLRP5, NLRP7, ZFP57 [32,33].

Discussion
Using GWMA, we found four patients with MLID (three BWS and one SRS). MLID was mostly revealed in LoM-type epimutations at the maternally methylated iDMRs, in accordance with previous studies [5,8,9,[34][35][36]. LoM predominance at the maternally methylated iDMRs could be explained by the fact that the majority of iDMRs are of maternal origin [37]. It is also possible that the maternally imprinted loci may be more susceptible to LoM than paternally imprinted loci in the establishment and maintenance of methylation during gametogenesis [34,38]. However, concurrent hypomethylation of both paternally-and maternally methylated iDMRs (MEG3 and PEG10, respectively) was observed in SRS1 in our study. Coexistence of paternal and maternal LoM was also reported in previous reports [8,39,40], which suggested the postzygotic origin of epigenetic defects, as opposed to occurrence in germ-line cells.
The involved iDMRs varied in our patients with MLID, supporting non-recurrent methylation defects in MLID [33]. All 12 additional iDMRs in our patients with BWS-MLID Table 3  were previously described in patients with BWS, SRS, and other IDs showing MLID [9,19,20,33,39,41]. Likewise, two additional iDMRs (MEG3 and PEG10) in our patient with SRS-MLID were also previously reported in patients with SRS showing MLID [16,42]. Among three patients with BWS-MLID, we found two hypomethylated iDMRs (PPIEL: Ex1-DMR, and FAM50B:TSS-DMR) in at least two patients. Aberrant DNA methylation at PPIEL and FAM50B have been associated with bipolar disorder, and intellectual disability, respectively [18,43]. However, the clinical consequence of these affected loci is uncertain because our patients with BWS-MLID did not exhibit these previously reported clinical features.

Mono-locus (N = 9) Multi-locus (N = 3) P-value
Our patients with BWS-and SRS-MLID did not show any other clinical features due to hypomethylation of other loci, implying an epi-dominant effect of one locus above the others in clinical phenotypes [8]. Likewise, most patients with MLID showed only clinical features related to the original IDs [8,34,44]. Nonetheless, several cases showed atypical phenotypes such as developmental delay and congenital anomalies [9,13,16,33,45]. This epigenotype-phenotype divergence in patients with MLID possibly reflects somatic mosaicism, that is, the degree of methylation disturbances at the critical CpGs within other iDMRs remained at a subclinical level in the target tissues [5]. However, the influence of additionally involved iDMRs cannot be completely excluded. As MLID-associated clinical signs may manifest as patients grow up, careful monitoring is warranted to determine the possible effect of additional disturbances. For example, patients with LoM at the PLAGL1 locus require monitoring for early-onset diabetes. Besides, routine screening for childhood cancers may be required in patients with LoM at the DIRAS3 iDMR, which is associated with imprinted tumor-associated genes [46].
In our study, no causative variants of MLID have been found in exome sequencing analysis in three probands with BWS who exhibited multiple LoM. Until now, variants affecting either zygotic or oocyte-derived trans-acting factors have rarely been reported in some patients with MLID or in their mothers [47]. In particular, recessive variants of ZFP57, which prevent demethylation of genomic iDMRs during early embryogenesis, have been identified in seven of 13 probands with TNDM showing MLID [13]. Maternal-effect variants affecting components of the subcortical maternal complex, which plays an important role in imprinting establishment, have also been reported in women with reproductive problems, such as hydatidiform mole and miscarriages, and in healthy women with offspring with MLID and IDs, including BWS [15,48,49]. Convincing evidence has been provided for the involvement of NLRP2, NLRP5, NLRP7, and KHDC3L in the etiology of IDs, while the role of other maternal-effect genes, including PADI6, OOEP, TLE6, UHRF1, and ZAR1 is not definitely established [48]. A detailed investigation in a large number of patients with MLID may provide further insights into the imprint acquisition and post-fertilization maintenance of imprinted DNA methylation.
This study has some limitations. First, we determined the methylation status in only blood leukocytes, and possible effect of tissue mosaicism could not be excluded. Second, our study focused only on the DNA methylation patterns analyzed in the 850K array. Methylation defects at other non-investigated regions may also involve the complex genome-wide methylation phenomenon. Third, due to the limited number of patients identified with MLID, and relatively short follow-up duration, the functional consequences of disruptions at other imprinted loci are still unclear. Nevertheless, this study was strengthened by the extensive analysis of the epigenome in deeply phenotyped patients with ID with proven epimutation. We have expanded the spectrum of epimutated iDMRs associated with MLID in patients with BWS and SRS.

Conclusions
Using GWMA, we identified MLID in three patients with BWS and one patient with SRS. Additional hypomethylated iDMRs varied in the number and degree of affected regions in individual patients. Patients with MLID did not show any clinically distinct characteristics. However, the long-term phenotypic impact of other disturbed imprinted loci remains to be elucidated via expansion of patient cohort.
Supporting information S1